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We have examined a class of cohesive-zone models of dynamic mode-I fracture, looking both at 
steady-state crack propagation and its stability against out-of-plane perturbations. Our work is an 
extension of that of Ching, Langer, and Nakanishi (CLN), who studied a non-dissipative version of 
this model and reported strong instability at all non-zero crack speeds. We have reformulated the 
CLN theory and have discovered, surprisingly, that their model is mathematically ill-posed. In an 
attempt to correct this difficulty and to construct models that might exhibit realistic behavior, we 
have extended the CLN analysis to include dissipative mechanisms within the cohesive zone. We 
have succeeded to some extent in finding mathematically well posed systems; and we even have 
found a class of models for which a transition from stability to instability may occur at a nonzero 
crack speed via a Hopf bifurcation at a finite wavelength of the applied perturbation. However, 
our general conclusion is that these cohesive-zone models are inherently unsatisfactory for use in 
dynamical studies. They are extremely difficult mathematically, and they seem to be highly sensitive 
to details that ought to be physically unimportant. 
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I. INTRODUCTION 

In a recent publication, Ching, Langer and Nakanishi 
(CLN) described a direct attempt to determine the 
linear stability of mode-I fracture against bending defor- 
mations. Their basic idea was to use a two-dimensional 
cohesive-zone model of the kind introduced by Dugdale 
lH and Barenblatt j|] , and to compute the change in the 
trajectory of a crack induced by an infinitesimally small, 
static, spatially oscillatory shear stress. In this way, they 
hoped to learn about the mechanisms that cause rough- 
ening of fracture surfaces and that limit the speeds of 
crack propagation |^,|| . The major part of CLN was de- 
voted to the mathematically difficult job of performing 
the fully elastodynamic calculation that is needed to de- 
termine this perturbed trajectory. CLN concluded that, 
for cohesive-zone models that do not include dissipative 
forces, the trajectories are strongly unstable. 

Our present investigation has focused, first, on a re- 
examination of the mathematical problems that emerged 
in CLN and, second, on an attempt to build into these 
models some dissipative mechanisms that might provide 
more realistic descriptions of fracture stability. Our re- 
sults have been disappointing. The fundamental assump- 
tion in CLN was that the Barenblatt-Dugdale cohesive- 
zone models, suitably generalized to apply to bending 
non-steady cracks, are robust, mathematically well-posed 
dynamical systems that incorporate much of the basic 
physics of fracture in solids. This assumption now seems 
to be incorrect. 

The cohesive-zone models that we have considered 
have the following properties. 

(1) The crack exists on a semi-infinite surface within an 
isotropic, linearly elastic solid. We shall discuss only two- 
dimensional cases in which the crack is a line that ends at 
the crack tip, and the surrounding elastic solid remains in 



a state of either plane stress or plane strain. The crack- 
opening displacements, that define the fracture surfaces, 
are the linear elastic displacements away from the crack 
line. Note that we exclude consideration of plastic de- 
formation either in the bulk or on the fracture surfaces, 
and that we also exclude nonlinear elasticity even near 
the crack tip. 

(2) Cohesive forces act between the fracture surfaces 
and are functions of the crack-opening displacements and 
their time derivatives. These forces vanish when the 
crack opening displacements become larger than some 
fixed range of the interactions; thus there should be a 
well defined cohesive zone near the tip of the crack within 
which these forces are nonzero. The cohesive forces may 
be nonlinear functions of the displacements and may con- 
tain dissipative, rate-dependent components; but all non- 
linearities and dissipative effects are confined within the 
cohesive zone. 

(3) Real materials cannot support infinite stresses; 
therefore the elastic stresses must be bounded at ev- 
ery point in the system. As in the original Barenblatt 
analysis, this finite-stress condition should determine the 
structure of the cohesive zone uniquely. It should be anal- 
ogous, for example, to the finite-stress boundary condi- 
tion at the free end of a moving string. 

Implicit in the statement of this third condition is the 
question of how literally to take the physics of the co- 
hesive zone. Our original idea was that we should not 
take this aspect of the model literally at all but, in- 
stead, should insist that the dimensions of the cohesive 
zone be the smallest relevant lengths in the theory, and 
that its dynamic response to external forces plays the 
role of phenomenological boundary conditions to be im- 
posed at the crack tip. In other words, our cohesive zone 
should be just a mathematical device for bridging the 
gap between macroscopic elastodynamics and atomic- 
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scale mechanisms at the crack tip. There are other possi- 
bihties, for example, in polymer fracture where we might 
develop a model of fibrils forming and breaking within 
an extended process zone in front of the crack. But, in 
such circumstances, we would want to add a great deal 
more detailed physics than we are prepared to consider 
here. 

We have discovered a number of interesting but not 
entirely reassuring properties of the cohesive-zone mod- 
els that satisfy conditions (1) through (3) above. Our 
dissipative models, in all cases where they turn out to be 
mathematically well posed, exhibit stable crack propaga- 
tion at small speeds and undergo a transition to insta- 
bility at larger speeds. We believe that the high-speed 
instability is associated with the tip-stress anomaly de- 
scribed by CLN. As pointed out in their paper, the tan- 
gential stress in all models satisfying the above conditions 
exceeds the normal (opening) stress everywhere along 
the crack trajectory, including at the tip, for all nonzero 
speeds of crack growth. This is a "relativistic" effect; 
that is, it is a result of the way the stress fields transform 
to a frame of reference that is moving with the crack tip 
at some fraction of the sound speed. This extreme form 
of the conventional Yoffe |^ argument seems, incorrectly 
as it turns out, to indicate that these models are man- 
ifestly unstable for any choice of the cohesive forces — 
that the forward direction is never preferred by the tip 
stresses. Our introduction of a dissipative term in the co- 
hesive shear stress does produce stability at small speeds. 
However, we have serious reservations about this result. 

One of our reservations is based on the fact that stabi- 
lization at small crack speeds in dissipative models, that 
turn out to be mathematically well-posed, depend on 
small-scale features within the cohesive zone, in violation 
of our intuition that such features should be irrelevant in 
physically sensible models. A specific, physically implau- 
sible aspect of our results is that increasing the strength 
of the dissipative term in the cohesive shear stress seems 
always to drive the response of the system in the direction 
of instability rather than stability. 

Our most serious concern, however, is the mathemati- 
cal fragility that we have found in these models. We have 
discovered that many, apparently credible choices of co- 
hesive forces lead to mathematically ill posed systems. 
The equations of motion often fail to have unique, phys- 
ically acceptable solutions. When acceptable solutions 
do exist, as mentioned above, they often depend sen- 
sitively upon apparently unphysical small-scale features 
within the cohesive zone. We conclude, therefore, that 
this particular class of cohesive-zone models — despite 
its popularity in recent decades — is not suitable for use 
in theories of dynamic fracture. At the very least, we 
believe that useful models will have to include plastic de- 
formation in extended regions outside the crack tip, and 
that they will have to contain enough degrees of freedom 
to describe the dynamics of blunting at the tip itself. 

This report is organized as follows. In Section ^ we 
summarize our reformulation of the main results of CLN. 



We have tried to make this Section readable by relegat- 
ing almost all mathematical details to a set of Appen- 
dices. Those Appendices include analysis which is new 
and possibly useful for further investigations, but wh ich is 
peripheral to the main thrust of this paper. Section 
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a report on our efforts to include dissipative mechanisms 
in the steady-state theory, an exercise that seemed to us 
to be a necessary part of our investigation for reasons 
that we explain there. 



Our principal conclusions emerge in Section IV where 
we examine the mathematical properties of the equation 
which determines how an advancing crack responds to 
bending perturbations. We first point out an unexpected 
mathematical failure of the CLN analysis. Next, we use 
a simplified mathematical version of the problem to show 
how adding a dissipative term to the cohesive shear stress 
may cure this problem. We then look in detail at sta- 
bility analyses for the several dissipative fracture mod- 
els whose steady-state properties are described in Sec- 
Here we point out that one perfectly plausible 
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model, with apparently all of the necessary ingredients, 
fails (like CLN) to produce a mathematically well de- 
fined response to bending perturbations. We conclude 
this Section, and the main body of this paper, by de- 
scribing another model that behaves in much the way we 
had expected from physical considerations. That is, it 
is stable at small crack-propagation speeds, becomes un- 
stable at larger speeds, and even, for some values of the 
dissipation parameters, undergoes a Hopf bifurcation so 
that the instability occurs at a well defined wavelength. 
Unfortunately, we have no reason to believe that this 
model is representative of any broad class of physically 
motivated fracture models, or that its behavior is a ro- 
bust feature of such such a class. 



II. SUMMARY OF PREVIOUS RESULTS 

We start this report by summarizing the strategy and 
initial results of CLN. In a number of cases, we present 
those results in forms that are different from the original 
versions. The reader should refer to CLN and the Appen- 
dices in the present paper for a more complete account 
of the calculations leading to these formulas. 



A. Geometry and boundary conditions 

The first step in any stability analysis must be the 
choice of a steady-state configuration whose stability is 
to be tested. For fracture problems, we need to choose 
a system that is translationally invariant along the di- 
rection of crack propagation. We then must decide upon 
boundary conditions. There are several ways to pose an 
elastodynamic boundary value problem so that its solu- 
tion resembles steady-state fracture. 
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We might, for example, imagine that a sohd containing 
a crack is truly infinite in extent and is loaded in such 
a way as to create a concentrated stress at the crack tip 
strong enough to break the material. If we hold the crack 
in place by some constraint and then release it, the sys- 
tem must eventually reach steady state in a finite region 
around the moving tip of the crack. This situation is not 
well-suited to a stability analysis, however, because the 
fracture surfaces continue to move away from each other 
indefinitely far behind the tip, and part of the elastic en- 
ergy initially stored in the solid is always being carried 
away to infinity. 

A second mathematical device, frequently used in the 
literature, is to apply moving tractions to the crack faces 
at some fixed distance behind the tip so that the crack 
advances steadily at the desired speed. The steady-state 
solution in this case does not have a wave that propagates 
to infinity; however this method of controlling the crack 
extension by adjusting the loading defeats our purpose 
of studying the stability of a freely advancing crack. 

Therefore, like CLN, we consider a crack moving from 
right to left along the centerline of an infinite elastic strip 
occupying the region (—00 < a; < +00, —W < y < +W) 
in the x,y plane. Far ahead of the crack, at a; — > —00, 
the strip is uniformly strained by an amount Exx = £Too 
(the strain tangential to the crack axis), Syy = snoc (the 
normal strain) and, for the moment, Sxy — 0. Rigidly 
clamped boundary conditions a.t y — ±W assure that in 
the steady-state all of the elastic energy initially stored 
in the strip must be dissipated at the crack tip. From 
the beginning of the analysis, we assume that the half- 
width W is very much larger than any other length scale 
in the problem (the wavelength of the perturbing stress 
or the length of the cohesive zone), thus we carry out 
most of our calculations in the limit W — *■ 00. However, 
there are several places where we need to reintroduce the 
length W. For example, the fully relaxed width of the 
crack must scale like W, and the stress-intensity factor 
that characterizes the forces transmitted to the crack tip 
is proportional to Vw. 

To avoid the complication of the fixed grip boundary 
condition, CLN used a mathematical technique devel- 
oped by Barber, Donley and Langer 0. They considered 
two different boundary-value problems whose solutions 
are equivalent in a certain limit. The first is the problem 
described in the last paragraph in which the elastic fields 
satisfy simple wave equations in the bulk and fixed-grip 
boundary conditions at a distant boundary. The second 
is defined in an infinite space, but the fields satisfy a 
driven massive wave equation with a small mass. If the 
driving force is proportional to the desired rigid displace- 
ment and the mass is proportional to 1/W, the solutions 
of the two boundary value problems are equivalent in the 
limit of very large W. 



B. Basic structure of the stability analysis 

Suppose that the applied strains Sn 00 and etoo are 
such that the unperturbed crack moves at speed v in the 
negative x direction along the centerline of the strip. We 
choose our units of time so that v is measured in units 
of the transverse sound speed. To examine the stability 
of this crack, we compute its steady-state response to 
a small (i.e. first order) external force that produces a 
spatially oscillating shear stress along the a;-axis: 



^{ext) 
-'xy 



(x,0) 



(2.1) 



The dimensionless stress Tixy is measured in units of twice 
the shear modulus, and im is the amplitude of the pertur- 
bation whose wavenumber is m. In CLN, Yl'xy*'^ is defined 
in the entire x, y plane so that, in principle, it can be the 
result of tractions applied at the edges of the strip or of 
material irregularities near the centerline. However, only 
the value of si'^f*'' on the centerline is relevant for the 
first order calculation presented here. 

The goal of the calculation is to compute the perturbed 
centerline y = Ycen{x) of the resulting fracture to first 
order in £„, that is: 



Y,,,,{x) ^Y„,e'^'- = XY{ni,v) 



(2.2) 



(Throughout this paper, symbols with carets such as 
im and Ym denote Fourier amplitudes.) xy is a com- 
plex steady-state response coefficient that depends on 
the wavenumber m and the average crack propagation 
speed V. If xy diverges at some v and some real value 
of m, then we would conclude that the system undergoes 
a change in dynamic stability at that wavenumber and 
speed. More generally, poles of xy in the complex m- 
plane are equivalent to stability eigenvalues. According 
to (2.2), poles in the lower half m-plane correspond to 



stable modes, and changes in stability occur when poles 
cross the real m-axis. 



C. Crack-Opening Displacements, Stresses and 
Cohesive Forces 

Let the functions U^^\x) be the normal displace- 
ments, relative to the local orientation of the centerline 
y = Ycenix), of the "upper" and "lower" [— ] fracture 
surfaces. Similarly, let the u!^\x) be the corresponding 



tangential displacements. Then 

Un 



^ N ^ N 



(2.3a) 



is the crack-opening normal (mode I) displacement; 



Us 



1 r 
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(2.3b) 



is the crack-opening shear (mode II) displacement; and 
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Ut=- 



U, 



[+] 



(2.3c) 



is the average tangential displacement. (We shall not use 
Ut explicitly for present purposes, but it is needed in 
principle for a complete description of the crack surface.) 

In a similar way, we define the normal, shear and tan- 
gential components of the stress fields evaluated at the 
fracture surfaces: 



(-^i ^en) ~t- (^, ^en) 



'^six) = — ^ {x, Ycen) + ' (a;, l^en) 



and 



St (2;) = — Sji (a;, Ycen) + {x, Ycen) 



(2.4a) 



(2.4b) 



(2.4c) 



As before, the subscripts N, S, and T denote tensor com- 
ponents rotated into the local orientation of the centerline 

Ycen ■ 

The cohesive stress acting between the two fracture 
surfaces has two compo nents , Y^cN and Scs, defined in 
analogy to ( ^.4a ) and ( ^.4b ) respectively. We assume 
that these are functions of the local displacements and, 
when we include dissipative terms, the time derivatives 
of these displacements: 



'Ycn{x) = Y^cNpNix), Un{x)]; 



(2.5a) 



l^csix) = 1:cs[Un{x),Un{x), Us{x),Us{x)]. (2.5b) 

In principle, we could include higher time derivatives in 
these equations. Moreover, in an even more general the- 
ory, this local definition of the cohesive forces may need 
to be replaced by a nonlocal form in which these forces 
depend on the crack opening displacements everywhere 
in the plastic zone. We see no compelling physical rea- 
son for either of these possible extensions of the present 
formulation. 

In our linear response theory, the shear displacements 
Us must be infinitesimally small quantities which van- 
ish when the crack is moving steadily along a straight 
line. For reasons of symmetry, these quantities can en- 
ter the normal cohesive stress only at second order and 
must therefore be neglected. Similarly, the cohesive shear 
stress must be linear in Us and its time derivatives; and 
it should not depend directly on Us but, rather, on the 
shear angle Q defined by 



Usjx) 
Un{x) 



tane(a;) = 6(2;). 



(2.6) 



Without loss of generality, therefore, we can write (2.51;) 
in the form 



^cS 



cn[Un{x), Un{x)] (e(x) -f 77e(x)) , (2.7) 



where 77 is a dissipation coefficient and Q denotes the time 
derivative of 8. Ou r ch oice of unity for the coefficient of 
the first factor in (2.7) corresponds to the CLN central- 
force assumption. That is, apart from the rate-dependent 
factor rjQ, the cohesive stress behaves like a central force 
acting between opposite points on the fracture surfaces. 

A general feature of the cohesive stress in these mod- 
els is that both its normal and shear components must 
vanish when the crack is fully open, that is, when Uj^ix) 
exceeds the range of the cohesive interactions. We define 
the cohesive zone to be that region near the tip of the 
crack, < x < £, within which the cohesive forces are 
nonzero. The length of the cohesive zone, £, is a dynamic 
quantity that depends on the state of motion of the crack. 
It appears frequently throughout this analysis. 



D. Steps in the CLN Analysis 

The CLN analysis starts by transforming the equations 
of linear elasticity into a frame of reference moving in the 
negative x-direction at a speed such that the tip of the 
crack is always at x' = 0, that is, x = x' + xtip{t), where 



Xfip {i) 



(2.8) 



This transformation into a non-uniformly moving frame 
is essential because it allows us to deal non-perturbatively 
with the various mathematical singularities that occur at 
the crack tip. In principle, the new frame should be cho- 
sen so that the tip is also a,t y' — in the moving x' , y' 
plane. To first order in the perturbation i^m however, 
this extra part of the transformation turns out to be un- 
necessary. To simplify notation throughout the rest of 
this analysis, we simply drop the apostrophes after per- 
forming the transformation. 

The next step in the CLN analysis is to write down 
formal solutions of the equations of linear elasticity sep- 
arately in the two regions of the x, y plane above and 
below Ycen (2;), and then to evaluate the unknown coeffi- 
cients that occur in those solutions by imposing boundary 
conditions on the centerline. Ahead of the crack tip, the 
centerline is purely fictitious and the boundary conditions 
are simply statements that the stresses and displacements 
must be continuous there. Behind the tip, on the other 
hand, these boundary conditions are statements about 
tractions on the fracture surfaces. That is, the stresses 
at the fracture surface, Y,iq{x) and £5(2:), must be bal- 
anced by the cohesive stresses in the cohesive zone and 
must vanish where the crack is open behind that zone. 

It is also necessary, at this stage of the analysis, to 
separate the problem into parts which are zeroth order 
and first order in the perturbation £,„. CLN introduced 
an extra subscript, or 1, to indicate the order at which 
each function was being computed. That will not be 
necessary here. All symbols with subscripts N are zeroth 
order in £„, and all with subscripts S are first order. 
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The combination of elasticity and boundary conditions 
produces a set of Wiener-Hopf equations that can be 
solved for the unknown stresses ahead of the tip and the 
unknown crack-opening displacements behind it. At ze- 
roth order in Sm, the equation that relates Un{x) and 
Yim{x) is a restatement of the conventional, steady-state, 
cohesive-zone model. As in the original analysis of Baren- 
blatt, the condition that the normal stress Ejv(a;) remains 
finite at the crack tip is sufficient to produce a unique so- 
lution. At first order in im, the equation involving Us{x) 
and Y^si^) contains the information that is needed to 
find i;). Again, the condition that Ss(0) remains 

finite appears to produce a mathematically well defined 
equation for Us- As we shall see, the situation is not so 
simple. 



(2.13) 



Equation ( p.9[ ) must be supplemented by the Baren- 
blatt condition (non-divergent stress at a; = 0): 



Swoo — 





1/2 ^ 


2T:h{v)W_ 


Jq 



dx 



Ycn[Un{x)] (2.14) 



where, 



(2.15) 



is the normal stress infinitely far ahead of the crack de- 
termined by the applied strains Eatoo a-nd Etoo- 

We also need to know that the energy dissipated per 
unit length of crack extension is 



E. Zeroth Order Problem 



frac 



dxT,cN[UN{x)] 



dUN 
dx 



(2.16) 



The zeroth order, steady-state problem is most conve- 
niently cast in the form of a nonlinear, singular integral 
equation for the crack-opening displacement U n (x) . (For 
a more detailed derivation, see Appendix A.) This equa- 
tion is: 



dUN 



1 



dx 'Kh{v I jQ 
or, equivalently, 



dy 



X 1 



y x-y 



YcN[UN{y)]; (2.9) 



Un{x) 



1 



Trb{v) 



TTb{v) 

e 

dy 





dyZ{x,y) Y.cN[UNiy)] 

Vv + V 



Because the linearly elastic material outside the cohesive 
zone is energy-conserving, the flux Tfrac defined in this 
way accounts for the entire conversion of elastic energy to 
fracture energy in this system. If T,cN includes dissipative 
effects, Tfra c wi ll be a function of th e cra ck speed v. If 
we use first (|2!|) and then ( |2.14| ) in (|1|), we find 



r/rac(w) 




dyT.cN{x) 



1 



y x-y 



^cN{y) 



dx \ W 

YcN[X) = — S 
X K 



(2.17) 



In 



Vy- 



where the first version of the right-h and side serves as a 
definition of the kernel Z{x,y). Eqs. ( |2.9| ) and ( 2.10| ) are 
valid everywhere in the region < x <^ W . 

Throughout this paper, all singular integrals such as 
that appearing in (2.9) denote Cauchy principal values. 
For simplicity, we have suppressed the argument Un = 
vdU/dx in writing the function Scat. The quantity 



The last equality restates the fact that all the elastic 
ScAr[C/jv(2/)], (2.10) energy initially stored in the strip is dissipated at the 
crack tip. In the limit w — > 0, or in the case where the 
cohesive stress is purely non-dissipative, Satoo must be 
equal to the Griffith threshold stress Sq: 



Snoo — Sg — 



W 



1/2 



(2.18) 



b{v) 



f3iv' 



r(A/3t-/3o); 



(2.11) 



plays an important role in much of what follows. Here, 



where Tq denotes the non-dissipative (threshold) fracture 
energy. 

It will be important to recogn ize th at, outside the co- 
hesive zone where £ <^ x <^ W , ( p. 10 ) shows the conven- 
tional behavior of the crack opening displacement. 
Only the fir st te rm in Z{x, y) contributes in this limit. 
Thus, with (|t|): 



/3f^l--; A^^l-.^; 



Pi 



1-y; (2.12) 



Un{x) 



2Vi 



and the parameter k is the square of the ratio of the 
longitudinal to transverse sound speeds. b{v) vanishes at 
the Rayleigh speed vr, b{vR) = 0. 

The upper limit of integration in ( ^.9|) is the length of 
the cohesive zone defined in the paragraph following 
Eq. (2/7). Let 5 be the maximum range of the cohesive 
force. Then £ is determined by the relation 



7r&(w) Jo \/y 



dy 



—p ^cN{y) 



8W 



■KK.b{v) 



1/2 



^NoaVx- 

(2.19) 



Note that all of these results so far are completely inde- 
pendent of our specific choice of cohesive stresses. 

We shall find it convenient to use a scaled version of 
the above equations. Indeed, these scaling transforma- 
tions are needed in order to uncover the mathematical 
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structure of the problem. Let Sq be some characteristic 
stress, perhaps the yield stress at the crack tip. Then 
define 



(2.20) 



and 



J:cn[Un] = ^o^cNpNix)] = ^o^cNiCiO]- (2-21) 

Again, we suppress the explicit dependence of CTcAt on the 
derivative of ^, but must remember in later ap plications 
that it may be present. With this notation, (2.9) becomes 




7<^cNm% 



(2.22) 



and the Barenblatt condition (2.14) is 



-'Nc 



k£ 


1/2 ^ 


So / 

Jo 


27rb{v)W 



a^NiCm (2.23) 



The relation ( 2.13 ) that determines combined with 
(2.20), becomes 



£ = S 



'!rb{v) 



(2.24) 



It is natural, therefore, to measure £ and other lengths in 
the units of 6. We return to this steady-state problem in 



Section [I] 



F. First-Order Problem 

The CLN result for the response coefficient XY{m,v) 
has the form 



= — Xy ^ (to, w) = im Aeoo + 'Do{'m,v) + 'Di{m,v) 



where 



(2.25) 



(2.26) 



It is useful for mathematical purposes — almost essential, 
in fact — to restrict our attention to the limiting case 
m£ <^ 1. That is, we look only at perturbations whose 
wavelengths are very much larger than £, the length of 
the cohesive zone. As discussed in Section ^, we expect 
£ to be smaller than any other relevant length in the 
system, including 27r/m. Thus this assumption seems 
reasonable for both mathematical and physical reasons. 
We shall see, however, that it limits our ability to explore 
the analytic properties of xyirn, v). 

In the small-m limit, the function Volm^v) can be 
written in the form [see Eq. ( B29| )] 



(2.27) 



The quantity ^^{v) is a function only of w^, even for cases 
in which the cohesive stress contains dissipative terms. 
Vi^iv) has the limiting value 



^ / X / K — 1 



1/2 



(2.28) 



It is a positive, monotonically increasing function in < 
w < Wij, and it diverges as v approaches the Rayleigh 

speed vji. 

The first two terms on the right-hand side of ( 2.25| ) 
ma ke no ref eren ce to the cohesive shear stresses defined 
in (2.5t) or ( |2.7| ). CLN pointed out that these first terms 
can be obtained by neglecting the cohesive shear stresses 
but constraining the elastic stress S5 to vanish at the 
crack tip — in effect imposing a "if// = 0" constraint. 
This first part of the response function therefore consti- 
tutes a generalization to non-zero velocities of the Cot- 
terell and Rice (CR) 1^ quasi-static theory of crack ex- 
tension. The quantity —Aeoo is preci sely t he CR "T 
stress"; and the factor {—imWY^^ in (2.27), when in- 



cluded in (2.25), produces a weak (VF-dependent) insta- 
bility for negative values of Aeoq. 

All of the effects of the cohesive shear stress ar e con - 
tained in the t hird t erm on the right-hand side of ( ^.25 ). 
In analogy to (2.27), and again in the small-m limit, we 
can write this quantity in the form [see Eq. ( B32| )] 



T>i{m,v) ^ ~m^£{-tmW)^^^J:NooVi{v) [ ^S.sHO]- 

Jo Vi. 

(2.29) 

Here, 'Di{v), like 'Do{v), is a fu nction only of and has 
the same limiting value ( 2.2g ) at w = 0. It too is a 
monotonically increasing function of v which diverges, 
more strongly than 'Do{v), a,s v ^ vr. 

The dimensionless quantity Scs is proportional to the 
cohesive shear stress given in ( p.7[ ). It is defined by 



Scs[am='JcN[Cmrcs[am 



where 



T^cS 



[a(0] - ( 1 + im£'^ 



rjv da 

£7 + 



(2.30) 



(2.31) 



Here, we have transformed the linear form 0+?70 into the 
moving frame and have replaced the shear angle by a 
new function denoted a(^). The latter quantity contains 
much of the interesting structure of the prob lem. Apart 
from multiplicative constants [shown in Eq. (B46)], a(^) 
is the ratio of the shear angle O to the displacement of the 
perturbed centerline Ycen', that is, it tells us how much 
shear is induced by the bending of the crack within the 
cohesive zone. Note that we have kept a term imrjv in 
( p. 31 ); the dissipation length rjv is not necessarily small 
compared to £, and thus may be comparable to 27r/m. 
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The equation that determines a(^) is deduced from 
the Wiener-Hopf equation that relates Us to Eg, sup- 
plemented by the condition that S5 (x) be n qn-di verRent 
at 2: = 0. [See the analysis leading to Eq. (B34).] The 
result is a linear, inhomogeneous singular integral equa- 
tion: 



III. STEADY-STATE, UNPERTURBED 
SOLUTIONS 

CLN based their calculations entirely on the standard 
Dugdale model for which 



Pi Jo d^' 



In 



ScsHO]- (2.32) 



This equation, in the non-dissipative case, is equivalent 
to Eq. (5.23) in CLN. We discuss it in detail in Section 

[V of the present paper. 

To conclude this summary, we rewrite (2.25) in the 
form 

— (m, v) = imAeoc + im{—imWy^''^^Noo'^o{v)^{m, v) 



where 

$(m, v) = I + im£ 



T^a{v) Jo 



(2.33) 



(2.34) 



Note that Xy^(™j is proportional to im. This is a con- 
sequence of translational symmetry in the y-direction; it 
is the slope of the perturbed trajectory, that is, imYm, 
and not its magnitude that responds to the perturbation 
im- Mo re imp ortantly, the second term on the right-hand 
side of ( 2.34 ), the term that derives from the cohesive 
shear stress, is explicitly proportional to im£. We trace 
this factor to the symmetry requirement that the inclu- 
sion of the cohesive shear stress, that is an odd function 
of the shear angle 8, can only affect the curvature of Ycen 
and not its slope. 

At first gl ance, it would seem that we should drop the 
mi term in (2.34) to be consistent with our small-m ap- 
proximation. In that case, we would have just $ = 1 
and would again recover the Cotterell and Rice theory. 
The main result of CLN, however, was based on the ob- 
servation that a(^) may diverge at small values of m£ 
and w, and that the resulting behavior of $(m, v) in the 
complex m-plane may completely change the conclusions 
of CR. For similar reasons, we have retained the factor 
p-mii{^-(, ) ^Yie integral kernel in Eq. ( 2.32|) , because 



this factor can play a role at first order in m£. Other 
small-TO corrections to the kernel are of order mv£, and 
thus are doubly small near m = v = 0. 



Finally, we remark that the ratio 'Di{v)/'Do{v) in (2.34) 
is unity at w = and diverges as v — > u^,. The func- 
tions 2?i(f) and I?o(f) are defined in detail in Eqs.( B3C| ) 
and (B45). For present purposes, we need to note only 
that this factor enhances the effect of the cohesive shear 
stresses at large velocities. 



So for < UNix) < (5, 
for Un{x) > S. 



(3.1) 



Here, Eq is the yield stress that we introduced in ( 2.21 ), 
and 6 is the range of the cohesive force defined in (2.13). 
Although Dugdale interpreted Eg as the plastic flow 
stress, there is in fact no rate dependent dissipation in 
this model; all of the stored elastic energy must be con- 
verted into the fracture energy Tq = SYiq. As a result, the 
crack can advance at any velocity in the range < v < vr 
at precisely the Griffith threshold given in (2.18). At 



higher driving forces, only propagation at the Rayleigh 
speed vr can occur with a fraction of the stored elastic 
energy radiated to infinity. 

As pointed out by CLN, the lack of dissipation in the 
steady-state model is problematic. A linear stability the- 
ory is an analysis of the first-order response of a system to 
some change in the driving force, and such a calculation 
may not make sense if only one driving force is allowed. 
The CLN procedure seemed permissible because the only 
perturbations being considered were shear stresses that 
bent the crack but did not accelerate it. Once we dis- 
covered other difficulties in CLN, however, we realized 
that we would have to make sure that those problems 
were not related to this unrealistic aspect of the steady- 
state behavior. Accordingly, we have examined several 
dissipative steady-state models, and summarize our re- 
sults in the next few paragraphs. In short, we find that 
adding dissipation to the steady state makes some in- 
teresting differences in the stability theory but does not 
make qualitative changes in its mathematical structure. 

To start, we need to translate some formulas pertaining 
to the D ugda le mo del in to the scaled variables introduced 
in Eqs. ( ^.2C ) and (2.21). Because the function o'cAf[C(C)] 
is unity throughout th e coh esive zone < ^ < 1, we can 
immediately integrate ( 2.22| ) twice to find 



C(0 = 2VC-(l-01n(i^ 



(3.2) 



We th en compute the length of the cohesive zone i from 
(in, using C(l) = 2: 



T:b{v)S 

2En ' 



(3.3) 



which tells us that ^ is a w-dependent quantity that is of 
order J/Eq at small speeds and vanishes at the Rayleigh 
speed. 

The conceptually simplest way of adding dissipation 
to the Dugdale model was suggested by Glennie |^ , who 
modified (^^) by writing: 



7 



(3.4a) 

or, in terms of the scaled variables, 



Here, 



0, 



A=^^: ~S = 



TTb{v) ' 



m > s. 

Trb{v)S 



(3.4b) 



(3.5) 



Using ( ^.4b|) in ( |2.22| ), we obtain an equation that can 
be solved analytically by standard techniques The 
conditions that C(C) must vanish at ^ = and be positive 
for ^ > are sufficient to determine the solution uniquely. 
A compact form of the result has been given by Freund 



cos-d 



^i?/7r+l/2 
(Wj^ 



where 



tani^ = vrA = 



i_ds_ (1 



(3.6) 



(3.7) 



For the most part, the properties of this solution are 
what one might expect for a model of this kind. Near the 
Griffith threshold, (2.18), the crack speed v rises linearly 



from zero as a function of the stress above the threshold 
with a slope inversely proportional to the dissipation co- 
efficient A. For large values of A, (that is, for large Aw or 
for V close to w^f where b{v) becomes small), the length of 
the cohesive zone £ goes to a finite value and the energy 
flux becomes 



^fractiv) « -^To A(t;) 



^To^. (3.8) 



This formula illustrates the behavior of such models in 
the limit of vanishing dissipation. When A is very small, 
V must jump quickly from to values near vj^ for a very 
small increment in T fract above the Grifflth threshold. 

There is, however, one fatal deficiency of this model. 
For any nonzero value of Aw, the stress diverges at ^ = 1; 
and we cannot allow divergent stresses in cohesive-zone 
models. The divergence is clear in th e fig ures published 
by Glennie and Freund. To see it from ( |3.6[ ), note that the 
integral over s is well defined at ^ = 1 for 7 > but the 
prefactor diverges there. Thus d(/d^ and, accordingly, 
crcAr[C(C)] s-re unbounded at the back end of the cohesive 
zone. We therefore must modify this form of the cohesive 
stress before using it in a stability calculation. 

We have examined two qualitatively different forms of 
the cohesive stress, both of which cure the unphysical 



singularity at ^ = 1. The first was motivated by the idea 
that the linear, velocity-dependent stress shown in ( 3.4b| ) 
probably ought to saturate at arbitrarily high opening 
rates. Thus we write 



d^ 



(3.9) 



for < C(^) < S. The parameters Ai and A2 are defined 
in direct analogy with (^|^) , with the "bare" parameters 
now named Ai and A2. The ratio A1/A2 determines the 
highest possible dissipative stress whereas the difference 
Ai — A2 is the effective dissipation coefficient for small 
opening rates. 

Another feature of the Dugdale and Glennie cohesive- 
zone models that may be unphysical is the discontinuity 
in the stress at the back end of the zone, ^ = 1. A simple 
way to achieve continuity in the cohesive stress without 
introducing any new para mete rs is to replace Dugdale's 
square-law stress function (3T) by a triangular function 
that vanishes at ^ — 1. That is. 



1 - 



(3.10) 



for < CiO < S- 

Neither of the latter two forms of the cohesive stress 



produces an analytically solvable version of (2.22), thus 
we have resorted to numerical methods. In the process, 
we have discovered that numerical solution of singular 
integral equations is a difficult and subtle enterprise. Ap- 
pendix G contains a summary of our strategy and some 
remarks about the mathematical nature of these prob- 
lems. 

Our numerical analysis indicates th at bot h of the new 
models, defined by Eqs. (3^) and ( 3.1C ), have non- 
divergent, physically acceptable solutions. In the saturat- 
ing case, (3^), the length of the cohesive zone £ vanishes 
and Tfractiv) reaches a finite limit at vr. This means 
that steady-state fracture is possible only for a finite 
range of driving forces in excess of the Griffith thresh- 
old. For the continuous cohesive stress shown in (3.10), 
on the other hand, the length of the cohesive zone £ in- 
creases to some limiting value and Tfracti^) diverges like 
{vr — v)^^ near vr as in the Glennie model with the dis- 
continuous cohesive stress. We show some representative 
graphs of Tfract{v) in Figure gll. 



IV. FIRST-ORDER SOLUTIONS 

We turn now to the main effort o f this investigation, 
an analysis of the integral equation (2.32) which should 
determine the linear response of the crack to bending 
perturbations. To start, we expand the kernel in ( 2.32| ) 



8 



To 




FIG. 1. The theoretical fracture toughness computed 
for three different normal cohesive stresses. The Glennie 
cohesive stress (curve a) with AEq = 1 produces results 
similar to the model with the continuous cohesive stress 
and the same value of A (curve b). The model with the 
saturating cohesive stress (curve c) with AiEq — 1 and 
A2S0 = 0-1 leads to a finite fracture toughness at v = vji. 



to first order in m£ and integrate once over ^. The result 
is: 



Pi Jo 



(4.1) 



where 



(4.2) 



and 



M(C,e') = (e-e')ln 




(4.3) 



The function 5cs[a(0] defined in ( |2.30| ). 

Although they never wrote their equation in this form, 
CLN were s tudy ing the equivalent of ( [4.1D for the Dug- 
dale model (3^) and ?7 = 0, that is, for 5c5[a(^)] — a(^), 
when they concluded that this fully dissipationless sys- 
tem was manifestly unstable. They reached this conclu- 
sion by noticing that, in the Dugdale case as shown in 



(4.4) 



Therefore, if one sets m = i> = in ( |4.l| ), so that Pt/Pi = 
1, then a = constant is trivially a null right eigenfunction 
of the operator on the left-hand side (i.e. it is a non- 
trivial solution of the homogeneous equation). Thus they 



expec ted t hat a(^) and, accordingly, the second term in 
(f> in ( ^.34 ), must have a singularity near the origin in the 
TO-plane for s mall v. By equating coefficients of S,^^^ on 
both sides of (4T), they found 



(3/4) 



(1 — 1/k)v^ — imi 



(4.5) 



which implied a zero in <i>(m,u) (a pole in xy) at ml = 
i{l — 1/k)v'^. This unstable pole in the upper half m- 
plane seemed to be consistent with the tip-stress analysis 
mentioned in Section ^ of the present paper. 

What CLN failed to realize is that the left-hand side 



of (4.1) is a singular integral operator that has null right 
eigenfunctions, not just for m = v = 0, but for a con- 
tinuous range of non-zero values of m and v. It is a 
well known property of singular integral equations such 
as (4.1) [see ||T^] that, like inhomogeneous differential 



equations, they can have both homogeneous and partic- 
ular solutions. That is precisely what happens here. The 
CLN version of (4.1) has no unique solution but, instead, 
has a continuous, one-parameter family of solutions ob- 
tained from a particular solution by adding an arbitrary 
amount of the null (i.e. homogeneous) solution. We see 
no physical criterion for selecting among these solutions, 
thus we conclude that this version of the linear stability 
theory is mathematically ill-posed. 

M. Marder 1 12 1 has pointed out to us that the existence 
of the physically acceptable null eigenfunction of the op- 
erator on the left-hand side of Eq. (4.1) means that the 
steady state solution whose stability we are studying may 
be non-unique. The null eigenfunction implies, at least 
in a linear approximation, that the steady state solu- 
tion can contain symmetry-breaking, shear deformations 
in the absence of the perturbing shear stress. To com- 
pute the amplitudes of these deformations, and to to find 
out precisely the conditions under which they exist, we 
would have to carry out a fully nonlinear analysis of the 
zeroth-order problem including shear crack opening dis- 
placements Us{x) as well as the normal displacements 
Un{x). If this interpretation is indeed correct, then what 
we have identified as a mathematical pathology in the lin- 
ear response calculation, is an equally pathological lack 
of definition of the steady-state behavior of these models. 

To illustrate the mathematical difficulties, we consider 
the special case of (LI) with, as above, iScs[a(^)] = a(^) 
and m = 0. This equation has the form 

c * «(o ^ ao^io - (1 - ^) f\'z{^,e)a{a = e^', 

Jo 

(4.6) 



where 1/ = 1 - Pt/Pi = (1 - 1/k){v'^/2). In Fig. |l^, 
we show null right eigenfunctions of Ci,, i.e. solutions of 
Ci, * anuuiO — 0, for various values of v. These functions 
anuuiS,) are all perfectly well behaved throughout the in- 
terval < ^ < 1; there seems to be no reason to reject 
them on either mathematical or physical grounds. 
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FIG. 2. Well behaved right null eigcnfunctions of Ci, for 
three different values of v. 



The addition of dissipation changes this situation com- 
pletely by producing an unphysical singularity in the null 
eigenvectors and thus a selection mechanism. However, 
this does not turn out to be the systematic and reliable 
selection mechanism that we had hoped to find. As a 
first step in seeing what happens, we look at the same 
Dugdale steady-state model, again with m = 0, but now 
with a nonzero value of the dissipation constant rj. That 
is, we study 



where 



•qv da 
1 d4' 



(4.7) 



(4.8) 



For numerical purposes, we find it best to eliminate a(^) 
in favor of t(^) by writing 



= — / de exp 
Vv Jo 



rjv 



(4.9) 



Note that we have fixed a(0) = 0. Our argument for 
doing so is analogous to one that Langer and Nakanishi 
used in discussing a model with a similar dissipative 
term. For nonzero 77, we must avoid a discontinuous jump 
in the shear angle, and therefore a delta-function spike 
i n th e shear cohesive stress at the crack tip. For rjv ^ i, 
(4.9) produces a thin boundary layer within which a(^) 
rises rapidly from zero at ^ = 0. Here is the first — but 
unfortunately not the last — occasion in which we shall 
see a length scale much smaller than £ playing a sensitive 
role in our analysis. As mentioned in Section |, we do 
not believe that small-scale features of this kind can be 
taken seriously as physically significant ingredients of our 
theory; thus we discount results in this regime. Never- 
theless, we must find out what happens mathematically. 

We have computed the null right eigcnfunctions of 
say TnuuiC)^ show several of them normalized to 
unity at ^ = in Fig. IV. Most notably, addition of 



FIG. 3. Right null eigcnfunctions of for two values of 
speed V and fixed rj ^ 5 (in our units 77 is a length and 
is therefore measured in units of 5 as all other lengths 
are). We have normalized these functions to unity at 
^ = 0. Note that the coefficient of the singularity at 
^ = 1 changes sign. 



the r] term in (4.7) produces power-law singularities at 
^ = 1 of the form (1 - S)-^ . In Fig. 0, we show the 
exponent 7, determined by fitting the functional form of 
Tnuii to a power law in some region around ^ = 1, as a 
function of -qv / £ for v = 0.1. This singularity becomes 
very weak at small values of •qvl£ and therefore makes 
our numerical analysis difficult in that region. By using 
smaller and smaller fitting regions, and correspondingly 
smaller numerical grid spacings, we have estimated 7 for 
values of rjvl£ somewhat less than those shown in the 
Figure. We believe that 7 remains positive all the way 
down to 'qv/£ = 0, implying that all of these null solu- 
tions for the stress are divergent and therefore physically 
unacceptable. The function ^{"qv/ £) seems not to depend 
sensitively on v. 
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FIG. 4. The exponents of the power law divergence of 
the solutions at ^ = 1 extracted by fitting the functional 
form of the discrete solution in an interval of fixed width 
near ^ = 1. 

We also have computed certain par ticular solutions 
of the inhomogeneous equation (4.7), the so-called 
"minimum- norm" solutions defined in Appendix C. We 
denote these by the symbol Tmin{C)i ^"^^ show several of 
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them in Fig. [V with parameters corresponding to the 
null solutions shown in Fig. IV. These particular solu- 
tions also have power-law singularities at f = 1 with the 
same values of 7 as th ose shown in Fig. IV . The accept- 
able solutions of (4.7), therefore, have the form 



t(0 



i(0 + C{v,ri) 



(4.10) 



where the factor C{v,t]) must be chosen so that the sin- 
gularities in Tmin and Tnuii cancel at ^ = 1. 




FIG. 5. The minimum norm solutions of the bending re- 
sponse equation with the Dugdale normal cohesive stress. 
The values of the parameters are identical to those in 
Fig. IV. Note that the coefRcients of the singular part 



do not change sign. 

The resulting solutions for r(^) or, equivalently, a(^) 
look much like those described by CLN near v = 0. The 
reason for this behavior, as shown in Fig. IV, is that 
the coefficient of the divergent part of r„„n, i-e. r^ing = 
lim^^i(l — TnuiiiOj vanishes at some value of v. 
Therefore, C{v,ri) diverges at that value of w, produc- 
ing a divergence in a(^). At and near that value of w, 
the selected solution t(^) is dominated by the null solu- 
tion, which we know is very nearly a constant near w = 0. 
However, the fact that we recover something like the CLN 
solutions near threshold does not mean that we also re- 
cover the CLN instability. As we shall see, the results 
of the stability analysis are sensitively dependent on the 
behavior of these weak singularities in the solutions of 
the integral equation, and crude approximations do not 
seem to work. 

With the understanding that a nonzero dissipation co- 
efficient f] may (or may not) determine a unique phys- 
ically acceptable solution, we return now to (4.1) and 
report results of numerical solutions with nonzero m and 
various choices of the normal cohesive stress as described 
In general, once we have found an ac- 



in Section III 



ceptable set of solutions for a(^) for variou s values of m 
and V, we can evaluate Xy^{''^t''J) in (2.33) and look for 
poles in the complex m-plane. Because mW must be 
indefinitely large, we can neglect the first term on the 
right-hand side of ( 2.33| ) (the "T stress") and consider 
only <I>(m,ti) defined in (2.34). Remember that the ze- 
roes of $(m, v) then correspond to the poles of XY{m^ v), 
and that a change in stability occurs whenever a pole 
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FIG. 6. The coefficients Tsing = limj^i(l — £,)''t{^) of 
the singular parts of the null and the minimum norm 
solutions for for a fixed 77 — 26. 



crosses the real m axis. We emphasize again that, be- 
cause of the weak singularities in a(^) at ^ = 1, these 
calculations become numerically difficult when rjv/i be- 
comes appreciably smaller than unity. This numerically 
troublesome region is also where we are least confident 
about the physical basis for our cohesive-zone models. 

The simplest non-trivial case is th e Du gdale model for 
the normal cohesive stress, i.e. Eq. ( |3.l[ ), o'cAf[C(C)] = 1> 
but with nonzero rj in th e cohesive shear stress factor 
Tcs[a(^)] defined in ( 2.31 ). This is the case in which we 
thought the stabilizing effect of rj should become appar- 
ent, and perhaps it does. Indeed, we now find stability at 
small crack speeds. Recall from our analysis of ( 4/7 ) that 
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FIG. 7. The dependence of $(m, w) on —im£ for rj ~ 5 
and two values of the crack's speed, one above and the 
other below the critical speed Vc{rj). 

the physically acceptable solution for a(^) with m = di- 
verges at a non-zero value of the velocity, say Vc{ii). As 
a result, the behavior of $(m, v) as a function of —im£ is 
radically different for values of the crack speed above and 
below this critical speed. It is clear from Fig. |^that, for 
V > Vc, the function <^(7n,v) has a zero in the unstable 
half of the m-plane, and that the system changes from 
stable to unstable as v increases through Vc- We plot 
the position of this zero of $(m, w), i.e. iTfipQiQ^^ as a 
function of crack speed in Fig. The function fc(?7) is 
shown in Fig. IV. Note that Vc decreases with increasing 
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v/vr 

FIG. 8. Zeros of $(m, v) in the complex m-plane plotted 
as —iiripoiei vs. v/vr for Dugdale normal cohesive stress, 
ry = (5 is fixed. These zeroes are located on the imaginary- 
m axis. 



77, which is perhaps not surprising. The parameter 77 oc- 
curs in our equations only in the combination rjv/i, and 
this term should provide the dominant w-dependence at 
small speeds. At larger speeds, of course, the "relativis- 
tic" dependence of the other terms in ^{m,v) start 
playing important roles and presumably are responsible 
for the transition to instability. Because this stability 
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FIG. 9. The critical velocity Vc as a function of the shear 
dissipation coefficient 77 for the Dugdale normal cohesive 
stress. 

transition occurs at ttz = 0, the interesting behavior takes 
place at values of \m£\ that are small enough to be con- 
sistent with our approximations. Thus we are forced to 
the conclusion that increasing the dissipation strength 77 
decreases near-threshold stability in this model. How- 



ever, the ratio r]v/£ is small everywhere in Fig. IV, and 
actually decreases as 77 increases. Thus, this apparently 
unphysical result occurs only in the region where the dis- 
sipation length 77V is small compared to the size of the 
cohesive zone £. If we simply ignore this small-?; region 
and look only at situations where rjv > £, then we re- 
cover something like the CLN instability, i.e. the laige-v 
behavior shown in Fig. IV. But the model then gives 



us no sensible description of the transition from stabil- 
ity to instability at small velocities. The next case that 
we examined produced our biggest surprise and our most 




v/vr 

FIG. 10. The real and imaginary parts of the zeros 
—impoiel of $(?7x, in the complex 77T,-plane for the 
model with rate-saturating dissipation in the normal co- 
hesive stress. The values of the various dissipation coef- 
ficients are fixed at 77 = (5, AiEq = 1 and A2S0 = 0.05. 



convincing evidence that this class of models is deficient 
in i mpo rtant ways. The normal cohesive stress shown 
in (BTG), for which the non-dissipative factor vanishes 
continuously at the back edge of the cohesive zone, ap- 
pears to us to provide a perfectly reasonable model of 
forces acting at a crack tip. Howeve r, w hen we tried 
to use ( 3.1C| ) in the integral equation (4.1), we found a 
null eigenfunction of of the operator on the left-hand side 
that, so far as we can tell, has no unphysical singularity 
anywhere in the interva l < C ^ 1- [We set a(0) = 
as in (4.£) and solved (4.1) for Tcs (O-l We conclude, 
therefore, that — just as in the much simpler case stud- 
ied by CLN — this model is mathematically ill-posed for 
purposes of stability analysis. 

We do not claim to understand what it means that an 
apparently well-posed physical model of dynamic fracture 
has a mathematically undefined response to infinitesi- 
mally small bending perturbations. Nor do we know 
whether we have accidentally found just one unusually 
pathological example, or whether many such models be- 
have in this way — perhaps all models with cohesive 
stresses that decrease continuously to zero at large sepa- 
rations between the crack faces. In any case, this example 
diminishes our confidence in this approach to fracture 
dynamics. The one class of models for which we have 
found somewhat sensible behavior is that for which the 
normal cohesive forces s atur ate at high opening rates, 
i.e. the model shown in (3^). Just as in the Dugdale 
model discussed previously, divergences at the back end 
of the cohesive zone provide a selection criterion, thus 
the model seems to be mathematically well posed. We 
have found that adding dissipation to the normal cohesive 
stress changes some but not all aspects of the behavior 
of the poles of xy {m, v) while keeping them in the region 
of the complex m-plane where \m£\ is small enough to be 
consistent with our approximations. Some typical results 
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FIG. 12. Same as for Fig. fV| except that t] = 2S. The 
imaginary part of the —irripoie^ is not plotted since it 
vanishes identically for these values of the parameters. 



are shown in Figs. |^ |^ and IV, where we plot the real 
and imaginary parts of —impoiei as functions of v for two 
different values of i] — S and rj = 26, and two sets of val- 
ues of Ai^2 chosen in such a way as to keep their ratio 
constant. As before, we find stability at small v and a 
transition to instability at a critical velocity Wc('7, Ai, A2). 
The values of r]v/i at the transition are 0.651, 0.510, and 
0.513 respectively. Thus we again seem to be in the phys- 
ically implausible region in which the microscopic length 
£ controls the behavior. For the smaller value of 77, the 
transition is a Hopf bifurcation which determines a char- 
acteristic wavenumber, i.e. a nonzero real part of m, at 
which the initial instability takes place. Increasing dissi- 
pation in the normal cohesive stress increases the critical 
speed thus making the crack more stable. For the larger 
value of 77, the transition occurs at m = and, ominously, 
at a smaller value of Vc- Thus once again we have a situ- 
ation in which increasing 77, i.e. increasing the resistance 
to shear deformations, reduces rather than increases the 
threshold for instability. 




0.8 0.9 1 



In summary, our results are the following: 

The CLN model, which is a simple Dugdale cohesive- 
zone model with no dissipation in either the normal or 
shear stresses, is not mathematically well posed for pur- 
poses of studying stability of mode-I fracture against out- 
of-plane bending perturbations. 

Similarly, a more complicated and apparently realistic 
model that contains both normal and shear dissipation, 
and for which the normal stress vanishes continuously at 
the back end of the cohesive zone, is mathematically ill- 
posed. The stability analysis does not make sense for the 
same mathematical reasons that caused the CLN analysis 
to fail. 

Addition of a simple linear dissipative shear stress 
converts the Dugdale model into a mathematically well 
defined system that is stable at the Griffith threshold 
and undergoes an instability at a higher velocity. How- 
ever, this behavior depends sensitively on the mathemat- 
ical structure of the cohesive zone at unphysically small 
length scales, and therefore seems to depend on artificial 
features of the model. One symptom of the artificiality of 
the model is that the critical velocity at which the crack 
becomes unstable decreases with increasing dissipation 
strength. 

Our fourth, most complex, model, which contains both 
shear dissipation and a nonlinear dissipative contribution 
to the normal stress, is mathematically well posed and 
exhibits physically reasonable behavior at least in some 
regions of its parameter space. However, like the dissipa- 
tive Dugdale model, it also loses stability with increasing 
strength of the shear dissipation. We have no reason to 
believe that this model is representative of any broad 
class of physically motivated fracture models, or that its 
behavior is a robust feature of such a class. 
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FIG. 13. Same as for Fig. IV except that AiSq = 2 and 
A2S0 = 0.1 so that their ratio is kept constant. 
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APPENDIX A: DERIVATION OF THE 
STEADY-STATE EQUATIONS 

The Wiener-Hopf equations for the steady-state theory 
are most easily written in a Fourier representation such 
that, for example, 



-t-oo 



dk 



^Nix) = / ^I]^(fc)e^''% (Ala) 

/ + 00 
dxT,M{x)e-'''''. (Alb) 
-OO 

In this notation, the zeroth order equation [CLN(3.9)] is 



Ejv(fc) = sl+^(fc) 



= 2^Swoo<5(fc)-F(fc)C/^+^(fc) 



(A2) 



The superscripts in parentheses (±) indicate functions 
that have singularities only in the upper (or lower) half 
fc-planes because they are Fourier transforms of functions 
that are nonzero only on the positive (or negative) a;-axis. 
Thus, the first line of (A2) tells us in Fourier language 

that the stress Tin{x) consists of a part Ti^j^\x) that is 
nonzero only ahead of the crack, a; < 0, plus the nor- 
mal component of the cohesive stress, 'E,(.n{x), which is 
nonzero for x > within the cohesive zone. The sec- 
ond line of ( [A^ ) relates these stresses to Un{x), which is 
nonzero only behind the crack tip, x > 0. The quantities 
Satoo and K are defined in Section HE. 
The Wiener-Hopf kernel F{k) is 



F{k) 



\2WJ 



b^{v)k^ 



1/2 



biv) (a^ + fc2)i/2. 

(A3) 



This expression is a simple interpolation between exact 
results in the limits kW oo (the only case of direct 
interest here) and kW — > 0. The quantity 



2b{v)W 



(A4) 



is infinitesimally small in the limit of very large W. In 
addition to producing a correct value for i^(0), its pres- 
ence here defines the analytic properties of F{k) in the 
complex k- plane. 

To solve the Wiener- Hopf equation, (|A2|), we first fac- 
tor the kernel (A3): 



with 

F^+\k) 



F{k) ^ /'(+)(fc)F(-)(fc), 



b{v){a + iky^^; F(-)(fc) = (a- 



(A5) 

iky/^ 

(A6) 



In terms of these factors, the solutions for the displace- 
ment and stress fields are, respectively: 



ikul^\k) 



1 



Fi+){k) 



Fi-){0) 



-ikAi%\k) 



and 



ifcEjv'(fc) = -F(")(fc) 



^Noo 

Who) 



+ ikAiJik) 



where 



.(+) 

-•cN 



[k') 



2ni k' ~ k±ie F(^){k' 



(A7) 



(A8) 



(A9) 



A quick way to obtain the condition for nonsingular 
stress, i.e. the Barenblatt condition, is the following. In 
(|A^), the factor F^~Hk) ~ fc^/^ for large k produces 
the usual |a;|~^/^ singularity in the stress at small jxj. 
Both terms inside the square brackets in (A8) are con- 
stants at large k; thus the Barenblatt condition is simply 
the requirement that these constants cancel each other. 
Specifically, 



-'Noo 



^ ^ ' 27r F(-)(fc) 
dx 



K 


1/2 ^ 


_27rb{v)W 


Jo 



^cn[Un{x)]. (AlO) 



Th e seco nd form of this result is the same as that shown 
in ( 2.14| ). We have redundantly inserted the length of 



the cohesive zone, £, as the upper limit of integration as 
a reminder that l^cNiUNix)] vanishes for larger values of 

X. 

An advantage of this technique is that we can substi- 
tute ( |A10| ) into (A7) to obtain an expression for uj^\k) 



from which the singularity has been eliminated: 



ik U^j^^ = 



1 



dk' 



(fc') 



(All) 



F(+){k) J 2Tr k' - k + ie F(-)(k' 

It i s now safe and useful to compute the Fourier transform 
of ( All ) and, in doing so, to take the limit W — + oo, 
i.e. a in the factors F(±)(A:). The result is Eq. (U). 
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APPENDIX B: DERIVATION OF THE 
FIRST-ORDER EQUATIONS 



Our starting point is the first-order equation 
[CLN(5.6)] that relates the shear components of the dis- 
placement and stress: 



Es{k,m) = T.[1.\k,m) + '{k,m) = 2TTE„^S{k - m) 
- G{k, m) f7^+^ (k) + L{k, m) C/^+^ (fc - m) Ym, (Bl) 



In analogy with ( A2), ' (fc, to) is the Fourier transform 



of the shear stress in the unbroken region a; < 0, and 
^cs' ^^"^ Fourier transform of the cohesive shear 



stress. 



The various ingredients of (Bl) are: 

Em = £m + imAeooym.; Aeoo = ^iVc 



(B2) 



where 



N^+\k) 



dp (p{p) 



p-k 
'''^ dp ip{p) 
kt- p-k' 

p^\m{p) qt{p)\ 



tan(^(p) - 4^ N ' 

mp) 



and 



(B9) 
(BIO) 



qfik) = (3f {k - k,+ ) {k - fc,_); ki± - ± ze; 



(Bll) 



The function which plays the role of the Wiener-Hopf 
kernel, the analog of F{k) in (AS), is 



G{k,m) 



qtv'^{k — to)- 



{k\iqt-q^), (B3) 



where 



qf^k^^^-ik^mf, 

ql EE fc2 _^;2^fc-TO)^ 
,,2 



V 



[k — mY 



(B4) 



Finally, the function which determines the coupling be- 
tween the oscillating shear stress and the unperturbed 
crack is: 



L{k,m) = 



2 ikqQ 
qtv^ 
'2^ 



1 - 



k — m 



(B5) 



k — m 
2 



{^)\k-m\[f{fi^k + m) - p\0, 



Ik 



Despite appearances, G(fc, to) and L(A:, m) are finite at 
k — m and w = 0. 

We now write the Wiener-Hopf kernel (pBSf) in the form: 



(B6) 



G(fc,TO) = (?( + )(/£, TO)G(")(fc, to), 



where the superscripts (±) have their usual significance. 
CLN find these factors to be: 



G^+\k) = 



Pi h{v) i{k - kR+) 
Pt [e + i{k~kt+)V/ 



2 exp[iV(+)(fc)], 



(B7) 



and 



qtik) = Pt {k - kt+) [k - kt-); kt± = ±—- ± le; 



l±v 



mv 

kR+ = ± ■ — ± le. 



vr±v 



(B12) 



(B13) 



The quantities kt+, etc. locate the branch points 
of G{k, to) in the complex fc-plane. Note that they are 
all proportional to mv. The quantities kji± are wave 
numbers for Doppler shifted Rayleigh modes, and are 
also of order mv. 

In the limit k ^ mv, the factors G*^*^ look very similar 
to the F'-^'^ given in ([A^): 



G(+\k,m)^^^{e + zk)'^^ 
Pt 

G^~\k,m) w {e-ikfl'^. 



(B14) 



This will be a very useful approximation. 

Given the (not necessarily in the approximate 

forms (B14)), we can solve the Wiener-Hopf equation 
(Bl)) for the shear stress and the shear displacement. 
The results are: 



-G^^\k,m) 



i{k — TO)i]^ {k, to) 

Em. 



G( '>{m,m) 



i{k ~ m) ^ (fc, m)Ym 



+ i{k-' to) k[,g'{k,m) 



■ (B15) 



i(k — m)U^^^^ {k, to) 

Em 



G(+){k,m) 



G(- 



(to, toI 



i{k~ m) A^^ (fc, m)Y„, 



i{k — m) (fc, to) 



(B16) 
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Here, 



1 



dk' 



2Tii k' ~ k ± ie 



L{k',m)u'l^\k' -m) 
G<P>{k\m) 

(B17) 



With this notation, (B19) can be rewritten as the for- 
mally exact expression for the response coefficient shown 
in ( p^ ): 



and 



dk' 



■<+) 



{k',m) 



2m k'-k±ie Gi-){k',m)' 



(B18) 



The next step is the analo g of the derivation of the 
Barenblatt condition in ( |A10 ). Just as we used the con- 
dition that the zeroth-order normal stress be nonsingular 
at the tip to determine in (2.18) andf in, for example, 



( 2.24 ), so it appears now that we can choose Ym so that 
the shear stress is nonsingular. As we shall see, however, 
the situation is not so simple as it seemed in CLN. 
On the right-hand side of (B15), in the limit of large 



k, the factor &~\k,m) is proportional to fc"'^/^ and the 
quantity in square brackets goes to a constant. Therefore, 
without regularization, the shear stress would diverge like 
Accordingly, we regularize the stress by requir- 
ing that the large-A: limit of the quantity in square brack- 
ets be zero, thus apparently fixing the value of Y„i. The 
result is: 



r,(-)( , fdk L{k,m)ul;\k~m) ^y 

= —G^ '(m,m) f Y„, 

. f dk E^,tVfc,m) 

'(m,m) / ."--^ ^ — -. 



(B19) 



In principle, we can solve (B19) for y,„. To do so, we 
define the function A{x) by 



Usix) = Un{x) Q{x) = A{x) Un{x) Ym 



(B20) 



A{x) needs to be defined only inside the cohesive zone, 
< a; < £. It is also useful to define the function us{x): 



Us{x) EE usix)Yme' 



U^+\k,m) = us{k~m)Ym. 

(B21) 



As argued in Section O, the cohesive shear stress is nec- 
essarily linear i n M g(x) or, equivalently, in A{x); thus we 
can write [see (2.7)] 



l]i+)(fc,m) -Sof™ / dxa,N[UN{x)]T,s[A{x)]e''^''-'^^^, 
Jo 

(B22) 

where, in the moving frame of reference, 

dA 

TcsiAix)] = (1 + irmjv)A{x) + rjv —. (B23) 

dx 



= -Xy ^ (to, -y) — im Aeoo + I?o (m, u) -f I?i (to, t;) 



(B24) 



where 
Voim, v) 

and 

I?i(to, w) = 



(B25) 



-So 



dk & ^ (to, to) 
2^ G(-)(fc,TO) 



x/ dxacs[UNix)]T,s[Aix)]e-'^'-"'^^. (B26) 
Jo 



It is useful to write 'DQ{in, v) in the form: 



f>o{m,v) 



dx h(x) Un{x), 



where 



hix)= I *2^:!(!I^L(fc,TO)e-(^-™)^ 
27r G(-\k,m) 



(B27) 



(B28) 



The function L{k,m), defined in (B6), is a homo- 
geneous function of k and to of order 2; that is, 
L{k,m) — L{k/m,,l). Similarly, &^\k^ra) = 



^/m&^\k/m, 1). Thus, h{x) is a function only of mx. 
For m£ <C 1, h{x) is very slowly varying on the scale 
£; that is, it changes appreciably only for x of order 
1/to ^ £. As a result, the overwhelmingly largest con- 
tribution to the integral in (B27) comes from the region 
outside the cohesi ve zo ne and, so l ong as mW remains 
large, we can use (2.19) to rewrite (B25) in the form 



'Do{m, v) = im{—imWy^^T,Noo 'Doiv). 
Here, Vq is a function only of v [see CLN(6.13)] 



(B29) 



{—im)^/^ 



-,1/2 



Kb{v) 



dkG^-\m,m) L{k,m)e-'^° 
2^ G(-)(fe,TO) [e + i(fc-TO)]3/2 



(B30) 



(The convergence factor e~''^° is not strictly necessary 
here.) Its limiting val ue a s w is given in ( 2.28 ). The 
easiest way to deriv e (|B30 ) is, in effect, to step backwards 
and replace ^/x in ( 2.19 ) by 



(a;) ^/x 



40F 



-ikx 



dk 



(e - i/c)3/2 



(B31) 
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where 9{x) is the unit step function. 

There is an important new point here. Despite the 
appearance of Un{x) in the formula for 'Do{m,v), the 
latter quantity is completely independent of what model 
we use for the cohesive forces so long as we consider only 
perturbations whose wavelengths 27r/m are much longer 
than £. 



The formula ( B26 ) for Vi also can be simplified for the 
case m£ <C 1. If we reverse the order of integration in 
( p326| ) , s o tha t it appears as a real-space integral in anal- 
ogy to (B27), then the cohesive shear factor limits the 
integration to the interior of the cohesive zone. As a re- 
sult, the relevant values of k in the integrand in ( B2(^ ) are 
much larger tha n my , and we can use the approximation 
( ^1^ to write ( ^26|) in the form: 



Viim,v) 



G^-\m,m) So 

-7^ C^cs[Un{x)] TcS 



A{x)] e'' 



(B32) 



At this point in our analysis, the remaining unknown 
ingredient of xvi'm, v) is the function A{x) appearing in 
(B32). According to ( [B20D , A{x) describes the bending 
of the cohesive zone that is induced by the perturbing 
shear stress. It is the central dynamical variable in this 
theory. 



We must compute A{x) by solving (Bl(: ) for the shear- 
displacement function Us{x). Equation (B16) is an in- 
homogeneous linear integral equation that determines 



[/^+^(fc,TO), which appears explicitly on the left-hand 
side and implicitly on the right-hand side via the factor 



Ti'"^\k, to) in Aj,g''(fc, to). This equatio n is b est rewritten 



by using the regularization condition (B19) to eliminate 
the explicit Em, so that the quantity in square brackets 
on the right-hand side automatically decreases like A;~^ 
at large k. The resuh is [CLN(5.23)]: 

dk' k' — m 



i{k — to) u^^\k — m): 



1 



Gi+^k,m) 



2Tr k' -k 



G(-)(fc',TO) 



So / dxacs[UNix)]Tcs[A{x)]e 


(+)/ 



— i(k' —ra)x 



-L{k',m)U];^{k' ^m) 



and 



where 



h{x,y) 



H{x) 



dyh{x,y) 
dy 



(B36) 



i{k- 



2TTi GW{k,m) 



)^ f dk' L{k',m) e-'^'^^' -^)y 
2tt k-k' -ieG(-){k',m)' 
(B37) 



In deriving (B34), we have integrated once by parts in 
the second term on the left-hand side. This is legal be- 
cause the upper limit of integration, i, is still redundant; 
d'cN[UN{x)] is defined over the whole positive x-axis and 
still may be assumed to vanish as smoothly as is neces- 
sary for a; > £. There is no difficulty at a; = so long as 
the stress remains finite there. 

To make further progress we must develop approxi- 
mate expressions for H(x) and K(x,y) by exploiting the 
fact that the wavelength of the perturbing stress is much 
larger than the length of the coh esive zone. We begin 
by noticing that the expression (B35) for K{x, y) can 
be obtained from that for h{x,y) in (B37) by replac- 
ing L{k',m) by 1. Both expressions can be rewritten 
by defining 



-(±) 



{x, m) 



dk e±*'=^ 
2^G'(±)(fc,TO)' 



(B38) 



Note that T^^\x,m) = for a: < 0. With the use of 
these auxiliary functions we can now rewrite (B37) as 



dx" L{x", TO)r(-)(a;' + x" + y - x). (B39) 



Here L{x, to) is the inverse Fourier transform of L( k, to) . 
The expression for K{x, y) can be obtained from ( [B39| ) 
by substituting a delta-function 5{x") for L{x",m) 



(B33) 



It is useful — in deed, essential — to invert the Fourier 
transform in (B33) and study this equation in x- space. 
We find: 



H{x) ^ — [A{x)Un{x)]- 
dx 

f' \ d 
^0 / dy —K{x,y) 
Jo 



acs[UN{y)]Tcs[A{y)] (B34) 



Here, 

K{x,y) = 



dy 

dk e'C^-™):^ 
27ri G'(+) (fc, m)J 2tt k - k' - ie G(-) (fc', to) ' 

(B35) 



dk' 



^-i(k'-m)y 



y) = e~'"'(^-J^) / dx'T'-+\x')r^-\x' + y- x). 
Jo 

(B40) 

Approximate expressions for r(^-'(a;,TO) can be derived 
for mx ^ 1 by ignoring the small k structure of 
&^\k,ra). Notable exceptions are the contributions to 
p(±) from the poles at fc = which become large as v 
approaches vr even for small to. We obtain 

r(-)(:E,TO) = ^ - ^fcfl-e''^^-" dte'''^-'\ 

(B41) 

The analogo us exp ression for r(+) {x , to) can be obtained 
multiplying (B41) by I3t/f3ib{v) and replacing kn- by 
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—kii+. In all the results reported here, we have con- 
sidered only velocities that are not close to the Rayleigh 
speed that the second term in (B41) becomes important. 

With the assumption that mv£ <C 1, and the observa- 
tion that both X and y arc s mall of order £ in (B35), we 
can use the approximations ( B14 ) to find 



K{x,y) ^ 



nb{v)Pi 



In 



~im{x—y) 



(B42) 



For reasons discussed in Section |1IF 
factor e-™(^- 



we have left the 



in (B42). Remember that the ap- 



proximations (B14) are valid for small mv, that is, for 
small enough v at an y value of m£. The remaining m- 
dependence in (B42) plays some role near the v = 
threshold. 



As in the calculation of Vq in (2.27), the kernel h{x, y) 
is a function only of mx and my, and the integral in 
( B36D is determined accurately, for ■m£ ^ 1, by values 
of y outside the cohesive zone. Therefore, in deriving an 
approximate expression for h{x, y), we may use the small 
mx approximation only for F^"*"^ (x, to). Returning to the 
Fourier transforms, we obtain 



h{x,y) 



h{v)(ii 
dk 



-im{x — y) 



L{k, to) 
27r G(-)(fc,TO) 



dx' 



Jk(x-x' -y) 



(B43) 



Since the int egral in ( B36 ) is dominated by large values 
of y, we use (2. IS) to find 



Hix) 



3m2 (-iTO)i/2 S 



Nc 



2G(-)(to,to) V^&(w) 



Viiv) 



(B44) 



where 



U/2 



-,1/2 



K h{v) 



Pi 



dk & \m^ra) L{k,m)e'' 



-ikO 



2^ G(-)(fc,m) [e-f i(fc"TO)]i/2 



(B45) 



In this case, the convergence factor is necessary. Before 
reversing the order of integration over k and y, which 
we have done here, we must close the fc-contour in the 
lower half plane in ( B43| ) because we want y > x - 
the exponential there. Then, in ( ]B45| ), the factor e 



X m 

ikO 



enforces this closure rule. Unlike the situation in (B30), 
we would get an incorrect result here if wc closed the 
contour in th e up per ha lf pl ane. We have chosen the 
prefactors in ( B44 ) and ( B45 ) so that, as in the case of 
'Do{v), the f uncti on I?i(w) depends only on v and has the 
same limit (2.2S) as u 0. 

Fin ally, we ma ke the scaling transformations shown in 
( ^.20D and ( |2.2lD , and define a{x) by 



Aix) 



m^{-irmTW£f/^i:Nc 
&^'> (to, to) So 



■Vi{v)a{0. (B46) 



The resuh is that (B34) becomes (2.32) 



APPENDIX C: NUMERICAL METHODS 

The most challenging numerical problems that we have 
faced in this investigation occur in the inhomogeneous, 
l inea r, singular integral equations of the form ( ^.32| ) or 
(4.1), which are to be solved for a(^) or, equivalently, 
Tcslo-iO]- We also have encountered nonlinear singular 
integral equations in the steady-state problem described 
in Section |l|. Here, however, the numerical problems are 
less severe; there are no weak singularities that need to be 
resolved in order to select physically admissible solutions. 

Our basic strategy has been to replace the unknown 
function, say a(^), by a piecewise-constant approximant 
that takes on values in intervals [^i-i, where = 
and ^Af = 1- We have used regular as well as variable 
interval le ngth s. Then an integral equation of the form 
[compare (4.1)] 



/(e)«(0 



(CI) 



can be evaluated at points — Zi ^ to become 

a system of N linear algebraic equations for the N un- 
knowns ai. We take the Zi = (^i-i -I- Ci)/2 to be the 
centers of the intervals. Thus, 



N 



B, 



(C2) 



where 



Cy = f{z,)S,, + / dC Q{z^, = g{z,). 



(C3) 



In general, G is a non-self-adjoint, complex matrix. 

Because the matrix G in all of our applications repre- 
sents a singular integral operator, it is ill-conditioned in 
the sense that at least one of its eigenvalues vanishes in 
the continuum limit, N —^ oo . The system of equations 



(C2) has a homogeneous solution in that limit; that is, 
(CI) has a solution with g{£^) = 0. In general, particular 
solutions also exist; therefore we can construct families 
of solutions by adding arbitrary amounts of homogeneous 
solutions to particular solutions. Of course, as described 
in Section [V, not all of these solutions are physically ac- 
ceptable. In well-posed problems, we expect to be able 
to exclude all but one of them. 

For any large but finite N, C may be a well- 
conditioned, invertible matrix, but one of its eigenvalues 
is generally too small for numerical purposes. To deal 
with this difficulty, we perform a singular-value decom- 
position 



of Cij to write it in the form 



k 



U^kWkVl, 



(C4) 
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where the Wk are the analogs of the eigenvalues of a well- 
conditioned, i.e. non-singular matrix. (The term "singu- 
lar" , as in "singular matrix" or "singular- value decompo- 
sition" , refers to non-invertibility of a matrix, i.e. to the 
existence of null eigenvectors, and not specifically to the 
singularity of an integral operator in the sense of Muskhe- 
lishvili 10.) The matrices U and V are unitary. They 
are unique up to permutations of columns (and corre- 
sponding elements of W) and their columns are left and 
right eigenvectors of C respectively. 



We then use the singular value decomposition (C4) to 
compute the so-called "minimum-norm" solution: 

ar"^J2'v^:l^Ul^B,, (C5) 

where the symbol means that, in the sum over j, we 
omit the term for which lim^v^oo Wj =0. In a finite 
matrix problem, the vector a™'" is an exact particular 



solution of (C2) if the vector B does not lie in the null 
space of C, and this particular solution would have the 
smallest possible norm, J^i 1^1 P- our case, although B 
generally does have a component in the direction of the 
r igh t null eigenvector of C, the minimum-norm solution 
( |C5| ) is guaranteed to be well behaved for large A^, and it 



does in fact converge to a particular solution of ( |Cl| ) . We 
have checked this convergence by computing the residual 
Si I C'y 1™" — Bip. This residual was generally of 
the same order magnitude as the smallest eigenvalue. It 
vanished in the limit of large TV as well. 

Using this procedure for the various cases discussed in 
this paper, we have been able to generate all members 
of what have always been one-parameter families of so- 
lutions. When both the homogeneous solution a""" and 
the minimum norm solution a™™ have been found to be 
singular at ^ = 1, the non-singular linear combination 
has been obtained approximately, up to corrections of 
order of the mesh interval squared, by the formula 

a, = a""" + — ^a""- (C6) 

We used a regular mesh with N = 400 to compute 
the zeroes of $(m,w) in the complex m-plane that are 



presented in Figures IV, IV, IV and IV. The detailed 



study of the ^ = 1 singularity in the null and minimum 
norm solutions of the bending response integral equation 



4.1 was done using a mesh with N — 800. The spacing 
of the mesh decreased as a power law near ^ = 1 to help 
achieve greater accuracy in resolving the singularity. 
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